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Abstract 

Exploratory simulations of four-nucleon interactions are performed taking into 
account the dynamical effects of internal nucleon loops. The four-nucleon interac- 
tions in the isoscalar and isovector channels are described by Yukawa interactions 
with auxiliary scalar fields. The nucleon mass and the average field lengths of 



the scalar fields are determined as a function of nucleon hopping parameter and 
Yukawa coupling strengths. There are no problems with "exceptional configura- 
tions" at strong couplings which make quenched simulations unreliable. 

1 Introduction 

There are impressive recent developments in solving QCD for the lightest nuclei, using the 
numerical technique of Lattice QCD (LQCD) (for recent reviews see [HI2])- In the QCD- 
approach the nucleons are described as bound states of quarks and gluons and the nuclei 
are composed out of these composite objects. For nuclei heavier than, say, the helium an 
approach based on quarks and gluons would require a tremendous computing power beyond 
the exaflops-scale. However, in order to compute, for instance, binding energies of heavier 
nuclei it is also not clear whether one has to stick to the fundamental theory, because the 
relevant degrees of freedom are not quarks and gluons but nucleons, pions and other hadrons 
described in a chiral effective theory [3], Hj. 

The chiral effective theory can be taken as a starting point for applying non-perturbative 
lattice thechniques for nuclear physics. A lattice approach to nuclear physics was presented 



in a row of papers by the authors of [5]-[TT]. For an overview over the various available 
approaches and methods to nuclear physics using chiral effective theories see the recent 
review articles [12]- [H]. The approach taken in these publications differs fundamentally 
from the methods applied in LQCD simulations where the path integral in Euclidean space- 
time is simulated using Monte Carlo methods. 

Recently, the application of Monte Carlo methods of LQCD has also been tried in [15] 
(for an earlier attempt see also [IS]). In Ref. [T5] two-nucleon binding energies have been 
calculated in the so called quenched approximation to the pionless effective field theory with 
isoscalar and isovector auxiliary scalar fields. In this approximation the fermion determinant 
in the path integral is replaced by a constant, i.e. internal fermion loops are neclected and 
hence no dynamical effects of the fermions are taken into account. Apart from being only 
an approximation to the Euclidean path integral, the quenched approximation suffers from 
the so called exceptional configurations - well known also from lattice QCD [171 [TS] . Such 
field configurations have extremely small eigenvalues of the fermion matrix and hence give 
huge contributions to some expectation values and completely spoil the statistics. As is 
has been concluded in Ref. [T5], in the effective field theory for nuclear interactions the 
occurence of exceptional configurations make quenched simulations impossible for strong 
couplings (actually for bare couplings in the range max(|Co|, |Ci|) > 0.3). This phenomenon 
has also been observed earlier in quenched simulations of some simpler Yukawa-models in 
[H]-[22j. Since this problem does not appear in numerical simulations of Yukawa models 
with dynamical fermions [231 [21], it can be expected that it does also disappear in nuclear 
Yukawa models if dynamical nucleons are included in the simulation update. 

In this paper the results of first unquenched Euclidean Monte Carlo simulations are 
presented in the model studied previously in [15] with particular emphasis on the strong 
coupling region max(|C |, |Ci|) > 0.3. After defining the lattice action and shortly describing 
the numerical simulation algorithm in Section [21 the results are presented. Conclusions and 
the discussion of the outlook for future work are contained in Section [31 

2 Numerical simulations 

2.1 Lattice action and updating 

The four-nucleon interactions are described here, as in Ref. [15] by Yukawa-couplings of the 
nucleon field to auxiliary scalar fields with isospin (4>x) an d 1 {<f>ax ( a = 1>2,3)). The 
nucleon fields are described by a pair of Grassmann variables and are denoted by ip ax and 
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ip ax where a = 1, 2 is the isospin index. The lattice action is defined by 

S = Sn + Sna ■ (1 
where Sn is the Wilson fermion lattice action [25J with the nucleon hopping parameter k: 

( _ ^= ±4 

x [ ^=±1 

and the scalar part of the action contains the Yukawa-type couplings to the nucleon field: 

Sna = £ KVf + 0iMJ + Co 0f $M + G x 0« farM] . (3) 

X 

Here C and C*i are the bare Yukawa-couplings, r a , (a — 1, 2, 3) are Pauli-matrices for isospin 
and a summation over repeated indices a is understood. The four-nucleon interactions are 
obtained after integrating over the auxiliary fields according to 

/oo „ 
d(j> exp{-0 2 -CftW)} = exp{— (^) 2 } . (4) 
-oo 4 

In Ref. [32] block fields were introduced in order to separate the physical cut-off due to 
the extended nature of hadrons from the inherent lattice cut-off given by the inverse of the 
lattice spacing a. In the present paper the block fields are not introduced and there is a 
unique lattice cut-off. The (straightforward) introduction of the block fields is left for future 
work. 

For the creation of sequences of interacting scalar fields 4>^\ 4>^1 the two-step Polynomial 
Hybrid Monte Carlo TSPHMC [261 EZl EH] algorithm is used which is a suitably adapted 
version of the original HMC [29] and, in particular, PHMC [30J algorithms. 



2.2 Numerical simulation results 

The numerical simulations are performed on a 8 3 ■ 32 lattice. This kind of lattices can be 
simulated on PC's (in our case on the PC cluster of the DESY theory group). The number 
of trajectories in the simulation samples is at least 5000, in some cases up to 15000. The 
nucleon mass and the average lengths squared of the scalar fields 

*.^X>f4°>), *^(^E««>. (5) 

are determined as functions of the three bare parameters k, Co, C\. (N x denotes the number 
of lattice sites.) 
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The results are summarized in Tables [10 and illustrated in Figures HE In general, 
the nucleon mass is increasing for increasing Cq and, somewhat surprisingly, decreasing for 
increasing C\. The decrease for increasing hopping parameter is expected since (2/t) -1 is, 
apart from an additive shift due to renormalization, the bare nucleon mass. This is well 
displayed by Fig. |3] which shows that for some critical hopping parameter k = k(Cq, Ci) the 
nucleon mass becomes very small. (It is expected that in the infinite volume limit it becomes 
actually zero.) This is important if one wants to perform simulations in the small lattice 
spacing limit. 

It is an interesting question what kind of behaviour is observed if the nucleon mass is 
divided by the expectation value of the average length of the isoscalar field (•I'o) 1 ^ 2 - As 
shown by Figs. HH5J for small Cq this ratio is approximately constant, but for larger values of 
Cq there is a moderate increase. For fixed Cq and C\ as a function of the hopping parameter 
the average field length does not change much. As a consequance, the qualitative behaviour 
in Figs. [3] and [6] is similar. 

3 Conclusion and outlook 

An important outcome of the performed unquenched simulations is the absence of any ex- 
ceptional configurations with very small eigenvalues of the fermion matrix. No such config- 
urations occured even at the strongest bare couplings which are much higher than the ones 
reached in Ref. [15J. The next goal should be to determine the renormalized four- nucleon 
couplings or, for instance, to tune the values of the dimensionless quantities cl Mn and 
A/m^v where Oq denotes the scattering length in the 1 S'o channel and A is the deuterium 
binding energy. This goes beyond the scope of this paper but is planned for future work. 

From the theoretical point of view it is an interesting question whether the tuning of 
the renormalized couplings to their physical values is possible at all by changing the values 
of the bare couplings Cq and C\. The model studied in this paper is a special case of a 
renormalizable Higgs- Yukawa model where the kinetic terms of the scalar fields are omitted 
and the bare quartic self-couplings of the scalar fields are set to zero. These terms are 
all dynamically generated by the four-nucleon interactions. This means that the model 
investigated in this paper explores a subspace in the space of renormalized couplings of the 
complete renormalizable Higgs- Yukawa model. Higgs- Yukawa models in general are expected 
to have a trivial continuum limit with all renormalized couplings equal to zero. Therefore, 
in the complete Yukawa-model and also in the subspace explored in this paper there is a 
lattice spacing (i.e. cut-off) dependent upper limit on the renomalized couplings which tend 
to zero if the lattice spacing goes to zero. (Zero lattice spacing means that the nucleon mass 
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in lattice units arriN is zero.) 

The model studied in this paper is, of course, incomplete and does not properly describe 
the nucleon interactions. First of all the pion field with its Yukawa interaction to the nucleons 
has to be introduced. Since the pion is even lighter than the nucleon, its kinetic term has 
to be included together with the quartic self interaction. For a complete description of the 
nuclear interactions further non-renormalized couplings have to be introduced which are 
generated within the framework of the chiral effective theory [3J H], which served as a basis 
of previous investigations of nuclear physics, for instance, in [5]-[14|. 

A further step towards a precise description of nuclear physics with Monte Carlo methods 
is to control lattice artifacts introduced by the discrete lattice which break Lorentz symmetry. 
One possibility is to introduce higher dimensional (non-renormalizable) couplings to com- 
pensate for the lattice artifacts within the framework of improved lattice actions |31j-[33j. 
Of course, in this way a large number of couplings has to be introduced and properly tuned. 
Therefore, this approach becomes at some point rather cumbersome. Another possibility 
is to introduce block fields [15] which mimic the extended nature of hadrons (nucleons and 
pion) and introduce by their finite space-time extensions a natural cut-off in momentum 
space. In this way the physical cut-off can be separated from the lattice cut off and a small 
lattice spacing limit can be defined keeping the physical cut-off at its desired physical value. 
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Tables 



Table 1: Numerical simulation results on 8 3 • 32 lattice for k = 0.10 and 
C\ = 0.05,0.10. Statistical errors in last digits are given in parentheses. 



K 


C 


Ci 


($o) 


($i> 

\ J- / 




0.10 


0.01 


0.05 


0.501442(86) 


1.48565(15) 


0.68519(10) 


0.10 


0.05 


0.05 


0.53377(10) 


1.48538(16) 


0.70653(37) 


0.10 


0.10 


0.05 


0.62883(12) 


1.48634(16) 


0.76671(40) 


0.10 


0.15 


0.05 


0.77137(15) 


1.48714(16) 


0.85493(57) 


0.10 


0.20 


0.05 


0.94305(17) 


1.48868(15) 


0.96010(72) 


0.10 


0.25 


0.05 


1.12759(19) 


1.48929(15) 


1.06971(72) 


0.10 


0.30 


0.05 


1.31323(20) 


1.49048(16) 


1.1813(10) 


0.10 


0.35 


0.05 


1.49272(22) 


1.49126(15) 


1.2927(10) 


0.10 


0.40 


0.05 


1.66176(23) 


1.49255(17) 


1.39602(80) 


0.10 


0.50 


0.05 


1.96256(27) 


1.49353(17) 


1.5840(14) 


0.10 


0.01 


0.10 


0.501517(84) 


1.44292(14) 


0.65804(19) 


0.10 


0.05 


0.10 


0.534201(93) 


1.44379(14) 


0.67974(28) 


0.10 


0.10 


0.10 


0.630929(42) 


1.446002(51) 


0.74205(18) 


0.10 


0.15 


0.10 


0.77541(16) 


1.44939(16) 


0.83309(79) 


0.10 


0.20 


0.10 


0.949356(51) 


1.453549(45) 


0.94291(20) 


0.10 


0.25 


0.10 


1.13598(18) 


1.45796(14) 


1.05644(79) 


0.10 


0.30 


0.10 


1.32251(21) 


1.46198(16) 


1.17402(64) 


0.10 


0.35 


0.10 


1.50271(21) 


1.46593(14) 


1.28402(95) 


0.10 


0.40 


0.10 


1.67156(24) 


1.46914(17) 


1.3872(10) 
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Table 2: Numerical simulation results on 8 3 • 32 lattice for k = 0.10 and 
C\ = 0.20,0.30. Statistical errors in last digits are given in parentheses. 



K 


C 


Ci 








0.10 


0.01 


0.20 


0.501340(96) 


1.28479(13) 


0.5629(12) 


0.10 


0.05 


0.20 


0.536649(99) 


1.28722(16) 


0.58162(80) 


0.10 


0.10 


0.20 


0.639301(76) 


1.296000(78) 


0.65369(54) 


0.10 


0.15 


0.20 


0.79218(16) 


1.30891(13) 


0.75732(85) 


0.10 


0.20 


0.20 


0.97365(12) 


1.32375(10) 


0.87765(60) 


0.10 


0.25 


0.20 


1.16692(22) 


1.33980(14) 


1.00478(60) 


0.10 


0.30 


0.20 


1.35767(24) 


1.35509(15) 


1.12985(85) 


0.10 


0.35 


0.20 


1.53994(24) 


1.36975(15) 


1.24829(94) 


0.10 


0.40 


0.20 


1.70999(27) 


1.38275(17) 


1.35859(87) 


0.10 


0.01 


0.30 


0.501390(95) 


1.06681(11) 


0.4296(47) 


0.10 


0.05 


0.30 


0.53978(12) 


1.07226(11) 


0.4501(25) 


0.10 


0.10 


0.30 


0.65177(14) 


1.08780(11) 


0.5431(36) 


0.10 


0.15 


0.30 


0.81674(18) 


1.11161(12) 


0.6495(19) 


0.10 


0.20 


0.30 


1.01076(22) 


1.14042(14) 


0.78928(71) 


0.10 


0.25 


0.30 


1.21305(26) 


1.17062(14) 


0.92980(98) 


0.10 


0.30 


0.30 


1.41081(29) 


1.20080(16) 


1.0690(10) 


0.10 


0.35 


0.30 


1.59643(29) 


1.22882(16) 


1.1968(13) 


0.10 


0.40 


0.30 


1.76870(31) 


1.25445(18) 


1.3087(37) 



9 



Table 3: Numerical simulation results on 8 3 • 32 lattice for k = 0.11 and 
C\ = 0.10. Statistical errors in last digits are given in parentheses. 



K 


Co 


Ci 


(<J>o) 


<*i> 




0.11 


0.01 


0.10 


0.501270(83) 


1.44445(14) 


0.4281(50) 


0.11 


0.05 


0.10 


0.533467(94) 


1.44488(14) 


0.4436(36) 


0.11 


0.10 


0.10 


0.62864(13) 


1.44667(15) 


0.5009(15) 


0.11 


0.15 


0.10 


0.77138(15) 


1.44977(15) 


0.6005(12) 


0.11 


0.20 


0.10 


0.94432(16) 


1.45419(15) 


0.72170(73) 


0.11 


0.25 


0.10 


1.13103(15) 


1.45808(13) 


0.85157(58) 


0.11 


0.30 


0.10 


1.31808(22) 


1.46208(16) 


0.97838(74) 


0.11 


0.35 


0.10 


1.49865(21) 


1.46600(15) 


1.10036(79) 


0.11 


0.40 


0.10 


1.66819(24) 


1.46932(16) 


1.21527(81) 
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Table 4: Numerical simulation results on 8 3 • 32 lattice for C\ = 0.05 and 
C = 0.40,0.60. Statistical errors in last digits are given in parentheses. 



K 


Co 


Ci 


L <*o> 


<<*>!> 


arriN 


0.10 


0.40 


0.05 


1.66176(23) 


1.49255(17) 


1.39402(80) 


0.11 


0.40 


0.05 


1.65810(23) 


1.49209(17) 


1.2230(13) 


0.12 


0.40 


0.05 


1.65354(24) 


1.49221(17) 


1.04973(95) 


0.13 


0.40 


0.05 


1.64713(23) 


1.49224(17) 


0.87697(84) 


0.14 


0.40 


0.05 


1.63763(23) 


1.49233(17) 


0.7033(13) 


0.15 


0.40 


0.05 


1.62421(25) 


1.49253(17) 


0.5154(20) 


0.16 


0.40 


0.05 


1.60449(25) 


1.49267(18) 


0.3612(94) 


0.10 


0.60 


0.05 


2.2115(24) 


1.49573(13) 


1.7512(36) 


0.11 


0.60 


0.05 


2.21459(29) 


1.49453(14) 


1.5983(31) 


0.12 


0.60 


0.05 


2.21269(31) 


1.49462(18) 


1.4567(11) 


0.13 


0.60 


0.05 


2.20967(33) 


1.49470(17) 


1.3159(11) 


0.14 


0.60 


0.05 


2.20585(32) 


1.49470(18) 


1.1725(19) 


0.15 


0.60 


0.05 


2.20051(32) 


1.49471(18) 


1.0368(14) 


0.16 


0.60 


0.05 


2.19391(27) 


1.49474(19) 


0.89612(85) 


0.17 


0.60 


0.05 


2.18459(29) 


1.49474(19) 


0.7542(11) 


0.18 


0.60 


0.05 


2.17155(28) 


1.49500(19) 


0.6109(21) 


0.19 


0.60 


0.05 


2.15482(30) 


1.49495(19) 


0.4664(40) 
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Figure 1: The nucleon mass as function of the coupling Cq for k = 0.10 
and different values of C\ . The lines are drawn to guide the eyes. 
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Figure 2: The nucleon mass as function of the coupling Cq for k = 0.10 
and k = 0.11. The lines are drawn to guide the eyes. 
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Figure 3: T/ie nucleon mass as function of the bare mass parameter (2k) _1 
/or C\ = 0.05 and /or Cq = 0.40 and Co = 0.60. The lines are drawn to 
guide the eyes. 
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Nucleon mass per scalar field length 
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Figure 4: The nucleon mass divided by the expectation value of the average 
length of the is o scalar field ($o) 1//2 a s function of the coupling Co for k = 
0.10 and different values of C\. The lines are drawn to guide the eyes. 
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Figure 5: The nucleon mass divided by the expectation value of the average 
length of the is o scalar field (^o) 1 ^ 2 a s function of the coupling Co for k = 
0.10 and k = 0.11. The lines are drawn to guide the eyes. 
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Figure 6: TTie nucleon mass divided by the expectation value of the average 
length of the is o scalar field (^o) 1 ^ 2 as function of the bare mass parameter 
(2k)- 1 for d = 0.05 and for C = 0.40 and C = 0.60. The lines are 
drawn to guide the eyes. 
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